set more off

use dataset, clear

keep if !mi(x)

*** do the graph and get the point estimate ***

DCdensity x, breakpoint(0) generate(Xj Yj r0 fhat se_fhat) nograph
local mytheta = `r(theta)'
local t = `r(theta)'/`r(se)'
su x
di 2*ttail(`r(N)', abs(`t'))

#delimit;

gr tw
	(sc Yj Xj, msym(O) col(gs10))
	(lpolyci Yj Xj if Xj < 0, ciplot(rline) alcol(black) clcol(black) alpat(dash) clwid(thick))
	(lpolyci Yj Xj if Xj > 0, ciplot(rline) alcol(black) clcol(black) alpat(dash) clwid(thick)) 
	
	,
		legend(off)
		plotregion(style(none))
		ylab(, angle(horiz))
		ytitle("Density")
		xtitle("Distance from state border")
		xline(0)
		xlab(-500(100)500)
		;
	

#delimit cr

gr export "_output/figureA5.pdf", replace
gr close

*** get the p-value: account for clustering by bootstrapping ***

gen B = .

su x

local myn = r(N)

forvalues i = 1(1)1000 {
	di `i'
	qui{
		preserve
			bsample, cl(state)
			DCdensity x, breakpoint(0) generate(Xj1 Yj1 r01 fhat1 se_fhat1) nograph
		restore
		replace B = r(theta) if _n == `i'
	}
}

su B

local myse = round(r(sd), .000001)
di 2*ttail(`myn', abs(`mytheta' / `myse'))
